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ABSTRACT 

Simulations predict that the first stars in a ACDM universe formed at redshifts z > 20 in minihalos 
with masses of about 10 6 M Q . We have studied their radiative feedback by simulating the propagation 
of ionization fronts (I-fronts) created by these first Population III stars (M* = 15 — 500M Q ) at 
z = 20, within the density field of a cosmological simulation of primordial star formation, outward 
thru the host minihalo and into the surrounding gas. A three-dimensional ray-tracing calculation 
tracks the I-front once the H II region evolves a "champagne flow" inside the minihalo, after the early 
D-type I-front detaches from the shock and runs ahead, becoming R-type. We take account of the 
hydrodynamical back-reaction by an approximate model of the central wind. We find that the escape 
fraction of ionizing radiation from the host halo increases with stellar mass, with 0.7 < / osc < 0.9 
for 80 < M*/Mq < 500. To quantify the ionizing efficiency of these stars as they begin cosmic 
reionization, we find that, for M* > 80M Q , the ratio of gas mass ionized to stellar mass is ~ 60, 000, 
roughly half the number of ionizing photons released per stellar baryon. Nearby minihalos are shown 
to trap the I-front, so their centers remain neutral. This is contrary to the recent suggestion that 
these stars would trigger formation of a second generation by fully ionizing neighboring minihalos, 
stimulating H2 formation in their cores. Finally, we discuss how the evacuation of gas from the host 
halo reduces the growth and luminosity of "miniquasars" that may form from black hole remnants of 
the first stars. 

Subject headings: cosmology: theory — galaxies: formation — intergalactic medium — H II regions 
— stars: formation 



1. INTRODUCTION 

The formation of the first stars marks the crucial tran- 
sition from an initially simple, homogeneous universe to 
a highly structured one at the end of the cosmic "dark 
ages" (e.g., Barkana & Loeb 2001; Bromm & Larson 
2004; Ciardi & Ferrara 2005). These so-called Popu- 
lation III (Pop III) stars are predicted to have formed 
in minihalos with virial temperatures T < 10 4 K at red- 
shifts z > 15 (e.g., Couchman & Rees 1986; Haiman, 
Thoul, k, Loeb 1996; Gnedin & Ostriker 1997; Tegmark 
et al. 1997; Yoshida et al. 2003). Numerical simulations 
are indicating that the first stars, forming in primordial 
minihalos, were predominantly very massive stars with 
typical masses M* > 100Af Q (e.g., Bromm, Coppi, & 
Larson 1999, 2002; Nakamura & Umemura 2001; Abel, 
Bryan, & Norman 2002). In this paper, we investigate 
the question: How did the radiation from the first stars 
ionize the surrounding medium, modifying the conditions 
for subsequent structure formation ? This radiative feed- 
back from the first stars could have played an important 
role in the reionization of the universe (e.g., Cen 2003; 
Ciardi, Ferrara, & White 2003; Wyithe & Loeb 2003; 
Sokasian et al. 2004). 

Recent observations of the large-angle polariza- 
tion anisotropy of the cosmic microwave background 
(CMB) with the Wilkinson Microwave Anisotropy Probe 
(WMAP; Spergel et al. 2003) imply a free electron 
Thomson scattering optical depth of r = 0.17, sug- 
gesting that the universe was substantially ionized by 
a redshift z = 17 (Kogut et al. 2003). Such an early 
episode of reionization may require a contribution from 
massive Pop III stars (e.g., Cen 2003; Wyithe & Loeb 



2003; Furlanetto & Loeb 2005). Analytical and numer- 
ical studies of reionization typically parametrize the ef- 
ficiency with which these stars reionize the universe in 
terms of quantities such as the escape fraction and frac- 
tion of baryons able to form stars (e.g., Haiman & Holder 
2003). In order to understand the role of such massive 
stars in reionization, it is therefore crucial to understand 
in detail the fate of the ionizing photons they produce, 
taking proper account of the structure within the host 
halo. 

Until now, studies of the propagation of the ionization 
front (I-front) within the host halo have been limited 
to analytical or one-dimensional numerical calculations 
(Kitayama et al. 2004; Whalen et al. 2004). These stud- 
ies suggest that the escape fraction is nearly unity for 
small halos, and as shown by Kitayama et al. (2004), 
is likely to be much smaller for larger halos. These con- 
clusions should not be taken too literally, however, since 
the structure of the halos in which the stars form is in- 
herently three-dimensional. Rather, that work should be 
viewed as laying the foundation for more detailed study 
in three dimensions. An effort along these lines was re- 
cently reported by O'Shea et al. (2005), focused on the 
dynamical consequences of the relic H II region left by 
the death of the first stars. 

Here, we present three-dimensional calculations of 
the propagation of an I-front through the host halo 
(M ~ 10 6 M Q ) and into the intergalactic medium (IGM). 
We model the hydrodynamic feedback that results from 
photoionization heating through the use of the similarity 
solutions developed by Shu et al. (2002; "Shu solution"), 
and calculate the propagation of the I-front by following 
its progress along individual rays that emanate from the 
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Fig. 1. — Hydrogen number density profile in a halo of mass 
M = 1O 6 M0 at z = 20. Solid line: Spherically-averaged density 
profile of minihalo in the SPH simulation. Dashed line: Density 
profile for a SIS with T SIS = 300 K. 

star. 

This paper is organized as follows. In §2 we present 
an analytical estimate of when and if the I-front escapes 
the host halo and describe our application of the self- 
similar Shu solution to primordial star forming halos. 
We describe our numerical method in §3, and present our 
results in §4. Finally, we discuss the implications of our 
calculations for reionization and further star formation 
in §5. 

2. PHYSICAL MODEL FOR TIME-DEPENDENT H II 
REGION 

At early times, as the I-front begins to propagate away 
from the star, its evolution is coupled to the hydrody- 
namics of the gas. The effect of this hydrodynamic re- 
sponse is to lower the density of gas as it expands in a 
wind, and eventually the I-front breaks away from the 
expanding hydrodynamic flow, racing ahead of it. In 
what follows, we describe an approximate, spherically- 
symmetric model for the relation between the I-front and 
the hydrodynamic flow in the center of the halo. This al- 
lows us to account for the consumption of ionizing pho- 
tons within the halo while at the same time tracking the 
three-dimensional evolution of the I-front after it breaks 
out from the center of the halo and propagates into the 
surrounding IGM. 

2.1. Early Evolution 

In a static density field, it takes of order a recombina- 
tion time for an I-front propagating away from a source 
that turns on instantaneously to slow to its "Stromgrcn 
radius" rs, at which point recombinations within bal- 
ance photons being emitted by the source. Generally, 
the I-front moves supersonically until it approaches the 
Stromgren radius, at which point it must become sub- 
sonic before it slows to zero velocity. The supersonic 
evolution of the I-front is generally referred to as "It- 
type" (rarefied), whereas the subsonic phase is referred 
to as "D-type" (dense). In the R-type phase, the I-front 



races ahead of the hydrodynamic response of the photo- 
heated gas. In the D-type phase, however, the gas is able 
to respond hydrodynamically, and a shock forms ahead 
of the I-front (see, e.g., Spitzer 1978). 

For the case considered here of a single, massive Pop III 
star forming in the center of a minihalo, this initial It- 
type phase when the star first begins to shine is likely 
to be very short lived, of order the recombination time 
in the star- forming cloud, t ICC < 1 yr for n ~ 10 6 cm -3 . 
This time is even shorter than the time it takes for the 
star to reach the main sequence, t ~ 10 5 yr, given by the 
Kelvin-Hclmholtz time. Thus, hydrodynamic effects are 
likely to be dominant at early times when the I-front is 
very near to the star. 

The study of the formation of these stars at very small 
scales defines the current frontier of our understanding, 
where the initial gas distribution and its interaction with 
the radiation emitted by the star is highly uncertain 
(e.g., Omukai & Palla 2001, 2003; Bromm & Locb 2004). 
For example, it is still not yet known whether a cen- 
trifugally supported disk will form (e.g., Tan & McKec 
2004) , or whether hydrodynamic processes can efficiently 
transport angular momentum outward, leaving a sub- 
Keplerian core (e.g., Abel et al. 2002). Omukai & Inut- 
suka (2002) studied the problem in spherical symmetry, 
making certain assumptions about the accretion flow and 
the size of a spherical H II region. Since the density and 
ionization structure in the immediate vicinity of accret- 
ing Pop III protostars is only poorly known, we must also 
make some assumptions here about the progress of the 
I-front at distances unresolved in the simulation we use 

(< 1 PC)- 

We will therefore assume in all that follows that the 
hydrodynamic response of the gas in the subsequent D- 
type phase will be to create a nearly uniform density, 
spherically-symmetric wind, bounded by a D-type I-front 
that is led by a shock. The degree to which the gas 
within this spherical shock is itself spherically-symmetric 
depends on the effectiveness with which the high interior 
pressure can reduce density inhomogeneities within. The 
timescale for this effect is the sound-crossing timescale, 
comparable to the expansion timescale of the weakly- 
supersonic D-type shock that moves at only a few times 
the sound speed. Since these two timescales are compa- 
rable, it is plausible that pressure will be able to homog- 
enize the density structure behind the shock. Further- 
more, the one-dimensional spherically-symmetric calcu- 
lations of Whalen et al. (2004) and Kitayama et al. 
(2004) indicate that pressure is capable of homogeniz- 
ing the density in radius, as shown by the nearly flat 
density profiles found behind the shock in the D-type 
phase. Although these are reasonable assumptions for 
these first calculations, we caution the reader that the 
detailed evolution of the I-front, especially very close to 
the star itself, can only be thoroughly understood by 
fully-coupled radiative transfer and hydrodynamic sim- 
ulations that resolve the accretion flow around the star, 
which we defer to future study. 
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Fig. 2. — Top: Density profile given by the Shu solution at 
different times after source turn-on, t = 5 X 10 4 ,10 5 , 2.5 X 10 5 , 
5 X 10 5 , 10 6 , and 2 X 10 6 yr, from left to right. Bottom: same as 
top but for velocity. The peak velocity, in the post-shock gas just 
behind the shock is constant in time, and is about 25% lower than 
the velocity of the shock itself. 

2.2. Model for Breakout 

Primordial stars are expected to form enshrouded in a 
highly concentrated distribution of gas. For a star form- 
ing within a halo with mass ~ 10 6 M Q , the spherically- 
averaged density profile of the gas, just prior to the on- 
set of protostellar collapse, is well-approximated by that 
of a singular isothermal sphere (SIS) with a temperature 
~ 300 K (Figure 1). For values relevant to a star- forming 
region in a mimhalo with mass M ~ 1O 6 M at z ~ 20, 
the hydrogen atom number density profile of the SIS is 
given by 



n(r) ~ 2.3 x 10 : 



( Jsis 
V300K J Vlpc 



(1) 



where we have assumed a hydrogen mass fraction X = 
0.75. We will take the SIS as a fiducial density profile for 
the calculations presented here. 

As discussed in §2.1, a D-type shock initially propa- 
gates outward just ahead of the I-front, leading to an 
outflow and corresponding drop in central density. After 
some time t = ie, however, the central density is suffi- 
ciently lowered so that recombinations can no longer trap 
the I-front behind the shock, and it quickly races ahead. 

This "breakout" time t B marks the moment at which 
ionizing radiation is no longer bottled up within and can 



escape. If the lifetime of the star t* < t B , then the front 
never escapes and the escape fraction is zero. This is 
essentially the reasoning used by Kitayama et al. (2004) 
to explain their result that ionizing radiation does not 
escape from halos with mass M > M cr it, where M cr ;t is 
determined by setting = t B . 

After breakout, the gas left behind is close to isother- 
mal with the high temperature of a photoionized gas (a 
few times 10 4 K). The strong density gradient results in a 
pressure imbalance that drives a wind outward, bounded 
by an isothermal shock. This "champagne" flow (e.g., 
Franco et al. 1990) has been analyzed through similar- 
ity methods by Shu et al. (2002), who found self-similar 
solutions for different power-law density stratifications 
p oc r~™, 3/2 < n < 3, and is also evident in the one- 
dimensional calculations of Whalen et al. (2004) and Ki- 
tayama et al. (2004). The family of solutions obtained 
by Shu et al. (2002) for the n = 2 case are described in 
terms of the similarity variables (eqs. (12) and (13) of 
Shu et al. 2002) 

x= — (2) 



and 



P(r,t) 



m p n(r) a(x) 



(3) 



X AirGt 2 ' 

where c s is the sound speed of the ionized gas and a(x) 
characterizes the shape of the density profile in the cham- 
pagne flow. If gas within the initial SIS has a sound speed 
c Sl i, then different solutions are obtained for a(x), de- 
pending on the ratio e = (c Si i/c s ) 2 , where c Si i and c s are 
the initial SIS sound speed and ionized gas sound speed, 



respectively. For Ti - 300K and T 
and the shock moves at v s = x s c s 



2xl0 4 if, e - 0.007 
t 40 km s _1 , where 
x s = 2.55. In Figure |3 we have plotted the density and 
velocity profiles in the Shu solution for the above param- 
eters. As seen in the figure, the density drops steadily 
in the center and is nearly uniform, while the velocity 
profile is unchanged as it moves outward. The peak ve- 
locity, corresponding to post-shock gas, is constant in 
time, ~ 30 km s -1 , and is less than the velocity of the 
shock itself. 

In what follows, we describe a model for when and 
where breakout occurs by finding the moment in the 
post-breakout Shu solution where the recombination rate 
inside of the shock is equal to the ionizing luminosity of 
the star. The condition for breakout is 

Q*=Aira B r 2 n 2 (r,t B )dr, (4) 

Jo 

where a B — 1.8 X 10 -13 cm 3 s _1 is the recombina- 
tion rate coefficient to excited states of hydrogen at 
T ~ 2 x 10 4 K, n(r,t) is the number density in the Shu 
solution, Q* is the ionizing photon luminosity of the star, 
and r B = r s h(t B ) = c s x s t B is the position of the shock 
at breakout. As shown by Bromm, Kudritzki, & Loeb 
(2001), the ionizing luminosity of primordial stars with 
masses M > 100M Q is roughly proportional to the mass 
of the star, Q* ~ 1.5 x 10 50 s _1 (M/100M Q ) (see Schaerer 
2002 for more detailed calculations). 

Combining equations J2Jl , © , an( i 0> , we can solve for 
the breakout radius 



Tb 



a B c s x s 



An(ji,mpG) 2 Q* 



l {x)x 2 dx. 



(5) 
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For fiducial values, we obtain 



rs " 2 - 3pc (ls) Gxi^-O • (6) 



Parameterizing the speed of the shock as v s , we can use 
the formula ts =TB/v a to derive the time after turn on 
at which breakout occurs, 



^5.6x1(^(40^) (3^) ( 3xl ^o s -i) 

Here and in the calculations we will present, we assume 
that the speed of the shock front in the D-type and cham- 
pagne phases is the same, v s ~ 40 km s _1 . The life- 
times of massive stars with masses 100 < M/M & < 500 
are within the range 2 < t* < 3 Myr (e.g., Bond, Ar- 
nett, & Carr 1984), much longer than our estimate of 
ts ~ 5.6 x 10 yr for our fiducial values. Thus, we ex- 
pect the time-dependent fraction of ionizing photons that 
escape from the halo to rapidly approach unity over this 
mass range. 

For lower stellar masses, and therefore lower values of 
Q*, the breakout time ts becomes comparable to the 
stellar lifetime i», which itself increases with decreasing 
mass (see Figure EJ. The precise value of the stellar mass 
at which t b — t* is therefore quite sensitive to the speed 
of the D-type shock, the density profile used in equa- 
tion and, of course, the main sequence lifetime of 
the star. In particular, the early hydrodynamic behavior 
of the gas in the D-type phase depends on an extrapo- 
lation to small scales where the mass distribution is not 
well understood. While our model for breakout is con- 
sistent with the one-dimensional calculations of Whalen 
et al. (2004) and Kitayama et al. (2004) in predicting 
that the escape fraction for massive stars M* > 100M Q 
approaches unity because breakout occurs early in their 
lifetimes, the threshold stellar mass at which tg = t* and 
the escape fraction goes to zero is not well determined 
and deserves further attention. 

3. NUMERICAL METHODOLOGY 

3.1. Cosmological SPH Simulation 

The basis for our radiative transfer calculations is a 
cosmological simulation of high-z structure formation 
that evolves both the dark matter and baryonic com- 
ponents, initialized according to the ACDM model at 
z = 100, to z = 20. We use the GADGET code 
that combines a tree, hierarchical gravity solver with 
the smoothed particle hydrodynamics (SPH) method 
(Springel, Yoshida, & White 2001). In carrying out the 
cosmological simulation used in this study, we adopt the 
same parameters as in earlier work (Bromm, Yoshida, & 
Hernquist 2003). Our periodic box size, however, is now 
L = 200/i -1 kpc comoving. Employing the same number 
of particles, Adm = Asph = 128 3 , as in Bromm et al. 
(2003), the SPH particle mass here is ~ 70M Q . 
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Fig. 3. — Timescales for breakout ts and stellar lifetimes t» 
versus stellar mass M*. Our estimate for tg becomes increasingly 
uncertain toward smaller stellar mass. Within these uncertainties, 
we estimate that no ionizing radiation will escape into the IGM for 
M* < 15M Q . 

We place the point source, representing the already fully 
formed Pop III star, at the location of the highest den- 
sity SPH particle in the simulation at z = 20, n ~ 10 4 
cm~ 3 , located within a halo of mass M v - lr ~ 10 6 A/ Q and 
virial radius r v j r ~ 150 pc. 

3.2. Ray Casting 

To calculate the propagation of the I- front, we con- 
struct a set of Nr = 12 x 2 12 = 49152 rays around the 
source. The directions of each of the rays are chosen to 
lie in the center of a HEALPixel 1 (Gorski et al. 2005). 
Each ray represents an equal solid angle of the sky as seen 
from the source. The rays are discretized into N s = 640 
segments of length Ar^ each, logarithmically spaced in 
radius, so that the radial coordinate of the front of the 
zth segment is r^. The inner radius is ~ 0.2 pc (proper), 
while the outer radius ~ 12 kpc (proper). In the case 
of a single source, we treat the evolution along each ray 
as independent of every other ray, so that we may repre- 
sent each ray as a different spherically symmetric density 
field. 

In order to calculate the evolution of the I-front, we 
must know the density of hydrogen along the ray. We do 
this by means of interpolation from a mesh upon which 
the density is precalculated. The density within each 
segment is assigned from the mesh at the midpoint of 
the ray segment, 

^+i/ 2 ^^(^+r 3 +1 ). (8) 

This discretization ensures that the midpoint of the ray 
is located at the point at which half the mass within 
the volume element is at r < f* i+1 / 2 and the other half 
is at r > r i+ i/2- The density value at the midpoint 
is determined by tri-linear intcrpolaton from the eight 

1 http:/ /www. eso.org/science/healpix 
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nearest nodes on the mesh, 

8 

n(x m , y m , Zm) = n g f(xg)f{y g )f{zg) 7 (9) 

where f(x g ) = 1— \x m — x g \/A c , A c is the mesh cell size, 
(x g ,y g ,z g ) are the coordinates of the eight nearest grid 
points, n g is the density at grid point g, 

x m =^+1/2 sin cos 4>, (10) 
Vm =r i+ i/ 2 sin0sin0, (11) 

~-m = n+l/2 COS0, (12) 

and (9, <p) are the angular coordinates of the ray. 

Given the substantial dynamic range necessary to re- 
solve the H II region around a Pop III star, the use of 
only one uniform mesh to interpolate between the SPH 
density field and the rays is not possible. Here we make 
use of the fact that the system is highly centrally con- 
centrated, which allows for the use of a set of concentric 
equal- resolution uniform meshes, each one half the lin- 
ear size of the last, centered on the star forming region 
in the center of the halo. Since the segments are spaced 
logarithmically in radius, the segment size at any point is 
smaller than the mesh cell of the highest resolution mesh 
that overlaps that point. We interpolate the SPH den- 
sity to each of the hierarchical meshes (see Appendix). 
For each ray segment midpoint, we find the highest res- 
olution mesh overlapping that point and use tri-linear 
interpolation from that mesh, as described above. 

3.3. Ionization Front Propagation 

In deriving the I-front evolution, we make the approx- 
imation that the front is sharp - i.e. gas is completely 
ionized inside and completely neutral outside. Because 
the equilibration time is short on the ionized side, ev- 
ery recombination is balanced by an absorption. Under 
this assumption, the I-front "jump condition" (Shapiro 
& Giroux 1987) implies a differential equation for the 
evolution of the I-front radius (Shapiro et al. 2005; Yu 
2005), 

dR = cQ(R,t) 

dt Q{R,t) + <±ivR 2 cn{Ry 1 ' 

where Q(R,t) is the ionizing photon luminosity at the 
surface of the front. This equation correctly takes into 
account the finite travel time of ionizing photons (e.g., 
R — » c as Q(R,t) — > oo). In general, this equation can 
be solved numerically, once Q(R, t) and n(R) are known. 

To approximate the hydrodynamic response due to 
photoheating, we combine the Shu solution with our ray 
tracing method to follow the I-front after it breaks out 
into the rest of the halo. We assume that the front makes 
an initial transition from R to D-type at radii r -C lpc, 
and creates a spherical D-type front that propagates out- 
ward to the breakout radius re, after which the Shu 
solution is expected to be valid. For each ray, we as- 
sume that the density profile of the gas at breakout is 
given by the self-similar solution inside of the shock, and 
is undisturbed outside, given by our cosmological SPH 
simulation (see §3.1). 

The initial I-front radius is independent of angle and is 
initialized to the breakout radius, so that equation 11311 
is solved with the initial value R(ts) = Tb- Q{R,t) in 
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Fig. 4. — Top: Instantaneous escape fraction for different masses, 
as labeled. Bottom: Time-averaged escape fraction (/ eS c) as de- 
fined in the text. Although the instantaneous escape fraction rises 
quickly just after breakout, the time-averaged value retains mem- 
ory of the breakout time and therefore lags behind. 

equation l|13fl depends on the density profile along a ray, 
and is given by 

Q(R,t) = Q* - 47ra B / n 2 {r,t)r 2 dr, (14) 
Jo 

where n(r, t) is given by the Shu solution at t for all 
r < r sh(t), while for r > r^t) the density is the initial 
unperturbed angle-dependent density distribution along 
each ray. We assume that the ionizing photon luminosity 
is constant over the lifetime of the star, with values given 
in Table 4 of Schaerer (2002). The final time in each run 
is set to the corresponding stellar lifetime for each mass, 
also given in Table 4 of Schaerer (2002) . 

4. RESULTS 

We have carried out several ray-tracing runs, each for a 
different stellar mass forming within the same host mini- 
halo. 

4.1. Escape Fraction 

The fraction of the ionizing photons emitted by the 
central star which escape into the IGM beyond the virial 
radius of the host minihalo is a fundamental ingredient 
in the theory of cosmic reionization and of the feedback 
of Pop III star formation on subsequent star and galaxy 
formation. We use our H II region calculations to 
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Fig. 5. — Mean escape fraction at the end of the star's lifetime t*, 
versus s tella r mass. Each symbol corresponds to (/ OS c) as defined 
by Eq. 11 (il for a different stellar mass calculation. Note that the 
escape fraction approaches zero for M < 15Mq, for which tg ^t t . 

derive this escape fraction / esc and its dependence on 
time during the lifetime of the star. Since our H II 
region density field and radiative transfer are three- 
dimensional, the escape fraction is angle-dependent. 
Along each ray, the escape fraction is given by 



/esc(^) — 



1 - 

0, 



J rvir n 2 (r,t)r 2 dr, R(t) > r v 
R(t) < r v 



( 15 ) 

where r V j r = 150 pc and n(r,t) is given by the Shu solu- 
tion for r < r s h(i), and by the SPH density field in that 
direction for r > r s h(t). The instantaneous escape frac- 
tion versus time, / CS c, determined by taking the average 
over all angles of the angle-dependent escape fraction is 
shown in the top panel of Figure 0] The average escape 
fraction between turn-on and time t < t„ is given by 



(fc 



f esc (t')dt', 



(16) 



and is shown in the bottom panel of Figure 0] Fig- 
ure El shows the average escape fraction at the end 
of the star's lifetime versus mass. For the very high 
mass Af» = 500Af Q case, the mean escape fraction is 
(/c SC ) ~ 0.9, while for M* = 80M Q it is about 0.7. We 
can understand the zero lifetime-averaged escape fraction 
at the smallest masses, as evident in Figure by com- 
paring the lifetime of the star and the breakout time. 
For ts < i* breakout occurs before the star dies and the 
escape fraction is expected to be greater than zero. For 
ts > <*, little or no radiation should escape. As can be 
seen in Figure E3 the threshold mass for which ts = i* 
is about 15 M . However, as we discussed in §2.2, the 
value of this threshold mass is very sensitive to the pa- 
rameters of our model. The escape fractions at masses 
M* < 50M Q are not robust predictions of our calcula- 
tions, but are shown here for completeness. 

4.2. Ionization History 



Shown in Figure is the evolution of the ionized 
gas mass outside the halo, Mnn(t), for different stel- 
lar masses. As expected, the more massive the star, the 
more gas is ionized. When expressed in units of the mass 
of the star, however, the quantity ?7hii = Mnu(t*)/M* is 
approximately constant with stellar mass for M* > 80, 
r/nii ^ 50, 000 - 60, 000 (Figure©. In the absence of re- 
combinations, so that every ionizing photon results in one 
ionized atom at the end of the star's lifetime, rjmi = Vph, 
where 

"* = "amT (17) 

is the number of ionizing photons produced per stel- 
lar H atom over the star's lifetime. This efficiency is 
roughly independent of mass for massive primordial stars 
M* > 5OM , ryph ~ 90, 000 — 100, 000 for X = 0.75 (e.g., 
Bromm et al. 2001; Venkatesan, Tumlinson, & Shull 
2003; Yoshida, Bromm, & Hernquist 2004). Recombi- 
nations cause the value of t/hii to be lower than rj p h by 
about a factor of two for large masses. 

4.3. IMF dependence 

Given the strong negative radiative feedback from H 2 
dissociating radiation that is expected once a star forms 
within a minihalo (e.g., Haiman, Rees, & Loeb 1997; 
Haiman, Abel & Rees 2000), it is unlikely that more 
than one star will exist there at any given time. This 
negative feedback may extend to nearby halos, though 
there is some uncertainty as to how strong this nega- 
tive feedback is (e.g., Ferrara 1998; Riccotti, Gnedin & 
Shull 2002). Thus, the first generation of stars forming 
within minihalos likely formed in isolation, and the initial 
mass function (IMF) of these stars was probably deter- 
mined by various properties of the host halos, such as 
their angular momentum and accretion rate. We make 
the reasonable assumption that the density structure of 
halos with a mass M ~ 10 6 M Q is universal, so that our 
determination of the escape fraction for this halo is close 
to what would be expected for other halos of compara- 
ble mass. For host halos of this mass, therefore, vari- 
ations in the escape fraction come only from variations 
in stellar mass. Under these assumptions, we can con- 
volve our results for one halo with different IMFs to see 
how the average escape fraction depends on the IMF. 
Usually, when applied to present-day star formation, the 
IMF describes the actual distribution of stellar masses in 
a cluster consisting of many members. In the primordial 
minihalo case, however, where stars are predicted to form 
in isolation, as single stars or at most as small multiples, 
the "IMF" would be more appropriately interpreted as 
a 'single-draw' probability distribution (Bromm & Lar- 
son 2004). Our analysis here is carried out in this latter 
sense. 

For definiteness, we use a Salpeter-like functional form, 
given by 



$(M) 



KM- 135 , 
0, 



M min < M < M n 
otherwise, 



and normalized so that 
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/ $(M)dlnM = l. 
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Fig. 6. — Mass ionized Mhii versus time for different stellar 
masses, as labeled. More massive stars produce more ionizing pho- 
tons in their lifetime and, therefore, ionize more of the surrounding 
gas. 

The total escape fraction, assuming one star forms per 
halo of mass M ~ 10 6 M Q , is given by 




IMF 

CSC 



f °° $(M)Q* (M)U (M)f csc (M)d In M 
f™ $(M)Q*(M)t*(M)dlnM 



(20) 



where the total number of photons released over a star's 
lifetime, Q^t* appears in the integrand because the es- 
cape fraction is being averaged over a period of time 
which is long compared to the lifetime of a star. For 



My, 



= 0.5 and M n 



500, which is a conservative 



estimate for the maximum Pop III stellar mass (Bromm 
& Loeb 2004), for example, the escape fraction is ~ 0.5, 
whereas for M m i n = 0.5 and M max = 80, it is ~ 0.35. 

4.4. Structure of H II region 

As can be seen from the visualization in Figure [8] 
the structure of the H II region is highly asymmetric, 
with deep shadows created by overdense gas. In partic- 
ular, nearby halos are not ionized, but rather are able 
to shield themselves and all that is behind them from 
the ionizing radiation of the star. This can clearly be 
seen in the bottom panels of Figure |H1 where overdense 
gas near to the central star remains neutral, despite be- 
ing so close. Figure shows the location of neutral and 
ionized SPH particles close to the star, showing that the 
highest density gas nearby the star remains neutral. For 
example, the highest density of hydrogen that is ionized 
within 500 pc of the 120M Q star is ~ 2 cm~ 3 at a ra- 
dius of ~ 200 pc, which corresponds to an overdensity 
of (5 ~ 4 x 10 3 , whereas the highest density of neutral 
hydrogen is ~ 400 cm~ 3 at a radius of ~ 220 pc, cor- 
responding to 6 ~ 2.5 x 10 5 . Similarly overdense gas 
that is further from the star is even more likely to shield 
itself and remain neutral, since the flux there is weaker 
because of spherical dilution. We find that ~ 4.9% of 
the sky at the end of the life of the 80M© star is cov- 
ered by high density gas that traps the I- front, whereas 
~ 2.6% of the sky is covered at the end of the 200A/© 



star's life. Such shielded regions are likely to be the sites 
of photoevaporation (Shapiro, Iliev & Raga 2004). The 
photoevaporation time of a 2 x 10 5 M© halo which is at 
a distance of 250 pc from a 120M Q star with luminosity 
1.4 x lO 5 ^- 1 is ~ 16 Myr (Iliev, Shapiro, & Raga 2005), 
longer than the lifetime of the star, so that most minihalo 
gas is likely to retain its original density structure. 

An important quantity associated with the "relic H II 
region" is the clumping factor. Shown in Figure lTUl is the 
clumping factor of the relic H II region, q = (n 2 )/n 2 , 
where the average is over the volume of the H II region 
and n is the cosmic mean density. Thus, the recombina- 
tion time in the H II region is given by t rcc = t r0 c,o/ c z, 
where i rcc ,o — 100 Myr is the recombination time of gas 
at the cosmic mean density. As the mass and luminos- 
ity of the star increase, the clumping factor of the relic 
H II region decreases. Apparently, clustering of mat- 
ter around the host halo causes the clumping factor to 
increase near the halo. Lower luminosity sources leave 
behind smaller, and therefore more clumpy, relic H II 
regions. The mean recombination time in the regions, 
however, is always less than the Hubble time ~ 175 Myr 
for even the largest stellar masses, with recombination 
times t rec < 60 Myr. These H II regions are thus likely 
to recombine unless other sources are able to keep them 
ionized. The timing of this recombination and the asso- 
ciated cooling of the recombining gas is crucial to under- 
standing the effect of photohcating on suppressing subse- 
quent halo formation, the so-called "entropy-floor" (Oh 
& Haiman 2003). 

4.5. I-front trapping by neighboring halos 

Whether or not nearby halos trap the I-front 
should determine whether ionization stimulates 
star formation in their centers. We can esti- 
mate the radius and density at which trapping 
occurs, as follows. The condition for trapping is 
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Fig. 7. — Ratio of ionized gas mass to stellar mass, »?hiIi versus 
stellar mass. Notice that for M* > 80AfQ, this ratio is almost 
independent of stellar mass. 
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Fig. 8. — Volume visualization at z = 20 of neutral density field (blue - low density, red - high density) and f-front (translucent white 
surface). Top row panels show a cubic volume ~ 13.6 kpc (proper) across, middle row ~ 6.8 kpc, and bottom row ~ 3.4 kpc. Left column 
is at the initial time, middle column shows simulation at t* = 3 Myr for the run with stellar mass M* = 80Mq, and the right column 
shows simulation at t* = 2.2 Myr for the run with stellar mass M, = 200Mq. The empty black region in the lower panels of middle and 
right columns indicates fully ionized gas around the source, and is fully revealed as the volume visualized shrinks to exclude the I-front 
that obscures this region in the larger volumes above. 



F 



a B nn(r)dr, 



(21) 



where F is the external flux of ionizing photons, r t is 
the radius at which the I-front is trapped, and r v i r is the 
virial radius of the halo. If we assume that the halo has 
a singular isothermal sphere density profile nn(r) cx r~ 2 



and an overdensity 5 

9(367r)V3Fm|02 



n 



, then solving for r t we obtain 

-, -1/3 



1 + 



a B X 2 M^ 3 (p(z)<5 vir ) 



s5/3 



"b J 



(22) 



where p(z) is the mean matter density of the universe 
at redshift z. For r^/r^ ir < 1, we can neglect the first 
term in the brackets, to see how this trapping condition 
depends upon the source and halo parameters and the 
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redshift, 



■>'i 



9(36 7 r) 1 /3 m 2 [ f72 i 



1/3 



(23) 

where po is the mean matter density at present. The 
density at the point where the I-front is trapped is 



Xp(z)n b S vil r 2 
n t = n H (r t ) = — 7z -■ 



(24) 



If we use the parameters of the minihalo nearest to our 
source halo for fiducial values, 



jO.18. 

r vir V 2 x 1O 5 M 



Q* 



V9 / r \ 2 /3 / r \ 5/9 



220 pc 



200 



1.4 x 10 50 s" 1 



-1/3 / 1 + Z \V3 
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(25) 



and 



nt ~ 3.6 cm 



2 x 10 5 M m 



-2/9 



220 pc 



-4/3 



(26) 



200 



-1/9 



1.4 x 10 50 s- 1 



2/3 



1 + Z 

21 



-1/3 



corresponding to a halo with total mass M v ; r = 2 x 
1O 5 M that is exposed to a flux F = Q»/(47rr 2 ) from 
a source with luminosity = 1.4 x 10 50 s _1 (for a 
stellar mass M* = 120Mq) at a distance r = 220 pc. 
Since M(< R) oc i? in a singular isothermal sphere, 
MHi/M V i r = r t /r V j r , and thus the neutral gas mass for 
the fiducial case above is ~ 5.5 x 10 3 M Q . The density at 
which the I-front is trapped is much smaller than the cen- 
tral density expected for a truncated isothermal sphere 
(Shapiro, Iliev, & Raga 1999), «h,o — 30 cm -3 at z = 20. 
Thus, nearby halos trap the I-front well before it reaches 
the central core. The weak dependence of r% /r v i r on halo 
mass and luminosity implies that trapping is a generic 
occurrence for halos surrounding single primordial stars. 



5. DISCUSSION 

We have studied the evolution of the H II region 
created by a massive Pop III star which forms in the 
current, standard ACDM universe in a minihalo of total 
mass M ~ 10 6 M@ at a redshift of z = 20. We have per- 
formed a three-dimensional ray-tracing calculation which 
tracks the position of the expanding I-front in every 
direction around the source in the pre-computed density 
field which results from a cosmological gas and N-body 
dynamics simulation based on the GADGET tree-SPH 
code. During the short lifetime (< few Myr) of such a 
star, the hydrodynamical back-reaction of the gas is rel- 
atively small as long as the front is a supersonic R-type, 
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Fig. 9. — Position ol selected SPH particles within 500 pc of the 
120 Mq star. Red particles are ionized and have a density above 
1 cm" 3 , all the neutral particles are shown in cyan, while only 
neutral particles with a density above 4 cm -3 are colored black. 
The radius of the shock in the Shu solution at the end of the star's 
lifetime, ~ 100 pc, is shown as the circle in the center. No SPH 
particles are shown within that radius. 

and, to first approximation, we are justified in treating 
the gas in this "static limit." At early times, however, 
when the I-front is still deep inside the minihalo which 
formed the star, the I-front is expected to make a 
transition from supersonic R-type to subsonic D-type, 
preceded by a shock, before it eventually accelerates to 
R-type again and detaches from the shock, racing ahead 
of it. 

To account for the impact of the expansion of the gas 
which results from this dynamical phase on the propa- 
gation of the I-front after it "breaks out," we have used 
the similarity solution of Shu et al. (2002) for cham- 
pagne flow. This solution allows us to determine when 
the transition from D-type to R-type and "break-out" 
occurs and, thereafter, to account for the consumption 
of ionizing photons in the expanding wind left behind 
in the central part of the minihalo. In this way, we are 
able to track the progress of the I-front inside the host 
minihalo and beyond, as it sweeps outward through the 
surrounding IGM and encounters other minihalos. This 
has allowed us to investigate the link, for the first time, 
between the formation of the first stars and the beginning 
of cosmic reionization on scales close to the stellar source 
that could not be resolved in previous three-dimensional 
studies of cosmic reionization. Among the results of this 
calculation are the following. 

Our simulations allow us to quantify the ionizing 
efficiency of the first-generation of Pop III stars in the 
ACDM universe as a function of stellar mass. The 
fraction of their ionizing radiation which escapes from 
their parent minihalo increases with stellar mass. For 
stars in the mass range 80 < M^/Mq < 500, we find 
0.7 < / osc < 0.9. This high escape fraction for high-mass 
stars is roughly consistent with the high escape fraction 
found for such high-mass stars by one-dimensional, 
spherical, hydrodynamical calculations (Whalen et al. 
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Fig. 10. — Top: Mean recombination time of H II region at end 
of star's life vs. stellar mass. Bottom: Clumping factor of H II 
regions at end of star's life vs. stellar mass. Less massive stars 
ionize a smaller volume, which implies a higher clumping factor 
because of clustering around the host halo. 

2004; Kitayama et al. 2004). For lower-mass stars, 
the escape fraction drops more rapidly with decreasing 
mass, as it takes a longer and longer fraction of the 
stellar lifetime for the I-front to end the D-type phase 
by reaching the "break-out" point, detaching from the 
shock and running ahead as a weak, R-type front to exit 
the halo. For M* < 15 - 20M Q , in fact, we find that 
the escape fraction should be zero and the I-front is 
D-type for the whole life of the star. More importantly, 
we find that this threshold mass is very sensitive to the 
hydrodynamic evolution of the I-front in the D-type 
phase. Given the great uncertainty regarding the inter- 
action of the stellar radiation and the gas immediately 
surrounding the star, a definitive answer to this question 
can only be obtained through three-dimensional gas 
dynamical simulations with radiative transfer that 
properly resolve the accretion flow onto the star. 

Once the H II region escapes the confines of the parent 
minihalo, the reionization of the universe begins. Our 
simulations yield the ratio of the final total mass ionized 
by each of these first-generation Pop III stars to their 
stellar mass. We find that, for M* > 8OM , this ratio 
is about 60,000, roughly half the number of H ionizing 
photons released per stellar atom during the lifetime of 
these stars, independent of stellar mass. 

We can obtain a very rough estimate of how effective 



Pop III stars are at reionizing the universe by assuming 
that all the stars have the same mass and form at the 
same redshift, each in its own minihalo of mass ~ 10 6 M Q . 
In the limit where the H II regions of individual stars 
are not overlapping, the ionized mass fraction fraction 
is given in terms of the halo mass function by /hii = 
Phii/ph where pun = rjunM^dn/dhiM is the mean den- 
sity of ionized gas, assuming each halo hosts a star of 
mass M*, and ionized a mass r]nu times the star's own 
mass, and pn is the mean mass density of hydrogen (For 
M* ~ 80M Q , for example, 7]u U ~ 50,000). Using the 
mass function of Shcth & Tormen (2002), dn/dhiM ~ 
130 Mpc~ 3 in comoving units at z = 20 for this mass 
range, we obtain / H n ~ 0.1[t/hii/5 x 1O 4 ][M*/8OM ]. If, 
instead, we use the ionized volume Vhii obtained directly 
from our calculations to derive the volume filling factor 
/v,hii = Vm\dn/dh\M for an 80M Q star, the final ion- 
ized volume is 7 x 10~ 4 comoving Mpc 3 , corresponding 
to a filling factor of /v,hii(M* = 8OM ) ~ 0.1. The sim- 
ilarity of the volume and mass ionized fraction indicates 
that the mean density in the ionized region is approxi- 
mately equal to the mean density of the universe. 

The above estimate is only the instantaneous ionized 
fraction, since the recombination times of each relic H II 
region are small fractions of the age of the universe at 
z = 20. This implies that many new generations of sim- 
ilar stars would have to form to continuously maintain 
this ionized fraction. A more conservative estimate of the 
effect of Pop III stars on reionization would also have to 
take account of the back-reaction of the starlight from 
earlier generations of stars on the star formation rate 
in the minihalos that follow (e.g., Mackey, Bromm, & 
Hernquist 2003; Yoshida et al. 2003; Furlanetto & Loeb 
2005). Since Pop III star formation depends upon the 
efficiency of H 2 formation and cooling inside minihalos, 
a background of UV starlight between 11.2 eV and 13.6 
cV is capable of suppressing this star formation by pho- 
todissociation the H2 following absorption in the Lyman- 
Werner bands. It is quite possible that the photodisso- 
ciating background builds up fast enough that minihalos 
are "sterilized" against further star formation before such 
a large fraction of the universe can be reionized in this 
way (Haiman et al. 1997). 

Hydrodynamic feedback due to photoionization heat- 
ing of the host halo will have a dramatic impact on its 
ultimate fate. Massive Pop III stars are expected to end 
their lives either by collapsing to black holes or exploding 
as pair-instability supernovae (e.g., Madau & Rees 2001; 
Heger et al. 2003). In both cases, it is important to 
know the properties of the host halo gas. Our model for 
the hydrodynamic feedback, in which an ionized, nearly- 
uniform density bubble bounded by an isothermal shock 
propagates outward at a few times the sound speed al- 
lows for an estimate of the density and size of the bubble 
at the end of the star's lifetime. For a lifetime of 2.5 Myr, 
the final size and density of the bubble are nibble — 100 
pc and ribubbic ^$ 1 cm~ 3 . If the star ended its life by ex- 
ploding, rather than collapsing to a black hole, then the 
SN remnant evolution inside this low-density bubble and 
beyond will differ from its evolution in the original undis- 
turbed minihalo gas. This should be taken into account 
in models of the impact of first-generation SN remnant 
blast-waves on their surroundings (e.g., Bromm, Yoshida, 
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& Hernquist 2003). 

This density can be used to estimate the accre- 
tion rate onto a possible remnant black hole, Mb — 
4 7 rG 2 M| H p/c 3 , (e.g., Bondi & Hoyle 1944; Springel, Di 
Matteo, & Hernquist 2005). Thus, e.g., for a black hole 
mass of IOOMq and a sound speed of 15 km s _1 , we ob- 
tain Ms ~ 2xlO _5 M Q Myr^ 1 . If we make the optimistic 
assumption that after a recombination time (~ 1.2 x 10 5 
yr) the gas can form H 2 molecules and cool back to ~ 300 
K without escaping from the halo, the accretion rate 
increases by a factor of 10 3 , to ~ 2 x 10" 2 M Q Myr _1 . 
These accretion rates are small compared to the Edding- 
ton accretion rate Mem — ^ttGAIbhttIp/ (earc) ~ 2M 
Myr" 1 , where the efficiency factor e = 0.1 (see also 
O'Shea et al. 2005). Such low accretion rates imply that 
remnant black holes from the first generation of stars 
are unlikely to be strong sources of radiation. Calcula- 
tions which do not explicitly take into account this re- 
duction in gas density near the remnant (e.g., Kuhlcn & 
Madau 2005) risk substantially overestimating their ra- 
diative feedback as miniquasars. These remnant black 
holes could begin to emit a substantial amount of ra- 
diation only after they encounter some other environ- 
ment containing cold, dense gas. It is not clear when or 
whether these black holes would ever encounter such an 
environment. At the very least, therefore, there should 
be some delay between the formation of the first genera- 
tion of stars and the X-ray emission from their remnants, 
if such emission ever occurs. 

Determining the fate of recombining gas in the relic 
H II regions left behind by the first stars is crucial. The 
contribution of these relic H II regions from the first- 
generation Pop III stars to the increasing ionized fraction 
of the universe during cosmic reionization depends upon 
their recombination time. Because of clustering around 
the host halo, the clumping factor and recombination 
time within the relic H II region depends on the mass 
of the star; higher stellar masses correspond to lower 
clumping factors and longer recombination times. The 
recombination time and clumping factor for a 4OM star 
are about 10 Myr and 12, respectively, whereas for a 
500Af Q star they are about 35 Myr and 3 (see Fig. 10). 

When ionized gas within the relic H II region cools 
radiatively and recombines, the nonequilibrium recom- 
bination lags the cooling, which enhances the residual 
ionized fraction at 10 4 K, promoting the formation of H2 
molecules which can cool the gas to T ~ 100 K and en- 
hance gravitational instability (Shapiro & Kang 1987). 
This corresponds to "positive feedback" if further star 
formation results (e.g., Ferrara 1998; Riccotti, Gncdin, 
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& Shull 2001). 

Recently, O'Shea et al. (2005) have addressed this is- 
sue of possible second generation star formation in the 
relic H II region of the first Pop III stars. They report 
that the I- front of the first star will fully ionize the neigh- 
boring minihalos and that, when the initial star dies, the 
dense cores of these ionized neighbor minihalos will be 
stimulated to form H2 molecules, leading to the second 
generation of Pop III stars. This assumes the initial star 
collapses to form a black hole without exploding as a su- 
pernova. We find, however, that the neighboring miniha- 
los are not fully ionized before the initial star dies, since 
the I-front is trapped by these minihalos and converted 
to D-type outside the core region, and the photoevapo- 
ration time for the minihalo exceeds the lifetime of the 
ionizing star. It remains to be seen, therefore, if the en- 
hanced H2 formation which O'Shea et al. (2005) found 
inside the nearest neighbor minihalo in their simulations 
of the relic H II region will occur in the presence of the I- 
front trapping and self-shielding which we report. More 
work will be required to resolve this question. 

Whether H2 molecules form in abundance or not de- 
pends on the density of the recombining gas. As we dis- 
cussed in §4.4, the highest density of gas in the relic 
H II region which we found to be fully ionized and, 
hence, capable of recombining to form H2 molecules, is a 
few cm~ 3 . A rough estimate of the H2 molecule forma- 
tion time is given by the recombination time of the gas, 
t lec ~ 10 5 yr. Even if H2 molecules form, however, it is 
not certain whether this will promote star formation in 
neighboring halos. As we have shown, the densest gas 
located in the center of nearby halos is not ionized. The 
gas that is ionized is not in the center, and is thus likely 
to recombine while being ejected from the halo as part 
of a supersonic, photoevaporative outflow (e.g., Shapiro, 
Iliev, & Raga 2004). Such gas is less likely to be grav- 
itationally unstable. In future work, we will investigate 
whether this gas is able to cool and form stars or simply 
gets evaporated into the diffuse IGM. 
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APPENDIX 

MASS-CONSERVING SPH INTERPOLATION ONTO A MESH 

Rather than interpolate directly from the SPH particles to our spherical grid of rays, we first interpolate the density 
to a uniform rectilinear mesh. The assignment of density to the mesh should conserve mass, which we accomplish as 
follows. We use a Gaussian kernel 



W{r, h) 



1 



-r 2 /h 2 



7T 3 /2/j3 

where h is the smoothing length and r is distance. This kernel is very similar to the commonly-used spline kernel, 



(Al) 







l>i. 



(A2) 



In our case, where a spline kernel has been used in the SPH calculation, we find that a Gaussian kernel is sufficient 
for interpolation purposes, provided we transform the smoothing lengths according to /iGauss = 7r _1 ^ 6 /ispiinc/2. 

For interpolation to a uniform rectilinear mesh with cell size A c , we wish to find the mean density within a cell 
centered at (x,y, z) contributed by a particle with smoothing length h located at the origin, 



W(r, h) 



1 



W(r, h), 



(A3) 



where the integral is over the volume of the cell. Since the kernel is a Gaussian, the spatial integral separates into 
three separate ones, so that 

= -Lz(x)E(y)E(z), (A4) 



where 



S(s) = erf 



8A3 
s + A c /2 



erf 



a - A c /2 



The value of the density averaged over a cell centered at r c is thus 

p(r c ) = ^2 m i W { r c ~*i,hi 



(A5) 



(A6) 



where the sum is over all particles i such that W(r c — r^, hi)jW(0, hi) > e, so as not to needlessly sum over particles 
with a negligible contribution. We find that a value e = 10~ 5 is sufficient for our purposes. 



